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THREE  DIMENSIONAL  TRANSIENT  NATURAL  CONVECTION 
IN  A  HORIZONTAL  CYLINDER:  A  NUMERICAL  ANALYSIS 


SUMMARY- 


This,  final  report  discusses  work  completed  under  USAF  Office  of 
Scientific  Research.  Contract  F^9620-79’’C-0l68  under  the  technical  guid¬ 
ance  of  Dr.  D.  G.  Samaras  in  the  period  16  May  1979  through  31  December 
1979. 

A  mathematical  formulation  of  the  governing  equations  for  transient 
natural  convection  in  a  finite  length  horizontal  cylinder  are  developed 
and  constructed  in  finite  difference  form.  The  boundary  conditions  con¬ 
sist  of  radial  heat  flux  for  a  specified  thermal  resistance,  axial  heat 
flux  from  one  closed  end  and  three  different  conditions  at  the  other  end 
to  represent  exposure  to  a  hot  convecting  gas  environment.  The  formula¬ 
tion  is  expressed  in  terms  of  the  vorticity  equations,  energy  equation 
and  a  set  of  vector  potential  equations.  Solution  is  by  the  ADI  (alter¬ 
nating  direction  implicite)  method  for  the  vorticity  and  energy  equations 
and  the  SOR  (successive  overrelation)  method  for  the  vector  potential 
equations. 

Numerical  experiments  were  run  using  the  model  to  determine  the  local 

wall  heat  flux  and  the  local  wall  temperatures.  Wail  thermal  resistance 

values  and  the  aspect  ratio  (length-to-diameter)  was  chosen  to  be  consis- 

tant  with  the  Air  Force  test  facility  at  Arnold  Air  Station.  A  heat 

transfer  correlation  is  presented  in  terms  of  the  Nusselt  and  Rayleigh 

numbers.  Steady  state  conditions  are  obtained  for  the  nondimensional  time 

(— y)  approximately  equal  to  .005.  Circumferential  heat  transfer  coeffi- 
o 

cient  variations  are  shown  with  larger  values  occurring  near  the  top  of  the 
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cylinder.  Axial  coefficients  vary  within  approximately  ±10%  with  the 
largest  values  occurring  near  the  center  of  the  cylinder.  With  respect 
to  test  conditions  at  the  Arnold  Air  Station  facility, the  convective 
components  appear  to  be  less  than  10%  of  the  radiative  heat  flux  to  the 
cylinder  walls  when  a  high  temperature  gas  (air)  is  enclosed  in  the 
cyl inder . 
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NOMENCLATURE 


A  -  representative  nondimensional  dependent  variable.  Equation  (.6) 

B  -  term  i,n  finite  difference  approximation,  Equation  (6) 

C  -  constant.  Equation  (6) 

0^  -  specific  heat  at  constant  pressure 

F  -  nondimensional  function  of  the  independent  variable 

g  -  acceleration  due  to  gravity,  directed  down 

G  -•  Grashoff  number  =  g8(2r  )3  T  /v2 

o  o 

G  -  Grashoff  number  =  GAT/T 


h 

i ,  j  ,  k- 
k 
L 
n 

Nu 
Pr 

<1 
O.i 

r 

r 

o 
R 

Ra 
S 
t 


convective  film  coefficient 

subscripts  in  the  finite  difference  approximations 
thermal  conductivity 
length  of  cyl inder 
time  step 


-  Nusselts  number  = 


h(2r) 


Prandtl  number  =  v/a 

heat  flux  at  a  solid  surface 

nondimensional  temperature  gradient  at  a  solid  surface 
radial  coordinate 
cyl inder  radius 

nondimensional  radial  coordinate  =  r/ro 
Rayleigh  number  =  Pr  Gr 

representative  nondimensional  dependent  variable,  Equation  (11) 
time  variable 


v 


<1 


T  -  temperature 

Tq  -  initial  temperature 

AT  -  prescribed  temperature  difference 

U  -  <(>  velocity  component 

V  -  r  velocity  component 

-  nond imens ional  velocity  vector 
w  -  z  velocity  component 

Z  -  axial  distance 

Greek  Letters 

6  ~  thermal  expansion  coefficient 

6  -  differential  operator 

0  -  nond imens ional  temperature 

v  -  viscosity 

t  -  nond imens ional  time 

$  -  circumferential  component 

iji  -  vector  potential 

to  -  nond imens ional  vorticity 


RESEARCH  OBJECTIVES 


The  specific  research  objectives,  of  this  study  are  as  follows: 

0)  Develop  a  model  of  three-dimensional  transient  natural  convection 
in  a  horizontal  cylinder  with  one  closed  end  and  one  open  end  which  is 
subjected  to  various  conditions  of  temperature  and  velocity. 

(.ii)  Perform  numerical  calculations  on  this  model  to  obtain  velocity 
and  temperature  distributions. 

(iii)  Integrate  the  foregoing  model  with  an  existing  radiation  model 
to  predict  the  inside  surface  temperature  of  the  cylinder. 
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STATUS  OF  RESEARCH  OBJECTIVES 


Mathematic  Model 

The  research  objectives  numerated  previously  have  been  completed.  An 
analytical  topi  for  examining  an  enclosed  transient  three--dimensional  cylin- 
crical  natural  convection  flow  field  has  been  developed.  The  impetus  for 
this  investigation  developed  from  the  need  to  predict  internal  temperature 
distributions  of  high  temperature  gas  containment  vessels  used  by  the  U.S. 
Air  Force  in  engine  test  facilities.  A  schematic  of  this  type  of  facility 
is  shown  in  Figure  1.  The  natural  convective  heat  flux  influence  on  the 
wall  temperature  history  appears  to  be  approximately  10%  of  the  radiative 
component.  The  test  facility  was  modelled  using  a  right  horizontal  cylin¬ 
der  lined  with  insulating  refractor  brick,  closed  at  one  end  and  open  at 
the  other  end.  The  coordinate  system  used  is  shown  in  Figure  2.  The  open 
end  is  in  direct  contact  with  hot  gases  and  has  been  modelled  under  various 
assumptions  discussed  later.  An  aspect  ratio  ( length-to-d iameter)  of  3-15 
was  used  in  ail  numerical  experiments. 

The  governing  equations  for  the  conservation  of  mass,  momentum  and 
energy  are  recast  in  terms  of  the  vorticity  transport  equation,  and  the 
energy  equation  with  a  defining  equation  for  the  vector  potential.  A 
similar  formulation  was  used  by  Aziz  and  Heliums  [1]  for  cartesian  coor¬ 
dinates.  In  vector  notation  the  system  of  equations  becomes: 


Propane  Fired  Heater 


Figure  2 


r 


Housing 

(hot  gases) 


Ambient  Conditions 


Geometry  and  coordinate  system  of  the  horizontal  extension 
arm  of  the  high  temperature  gas  heater. 
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(2) 


£»  -  i-  We 

Dt  Pr 


where  the  vorticity  is  defined  by: 

v  x  7  =  Z 


(3) 


such  that 

y2ii  =  -w  (*>) 

and  the  vector  potential  is  defined  by  the  relationship 

V  =  -V  x  *  (5) 

The  variables  used  in  the  above  equations  have  been  nond i mens ional i zed  with 
respect  to  the  time  scale  (r^/v) ,  and  length  scale  r^. 

This  mathematical  formulation  is  based  on  the  following  assumptions 
applied  to  the  fluid  system: 

(i)  Newtonian  behavior  with  constant  viscosity  and  thermal  conductivity 

(ii)  the  Boussinesq  approximation  where  the  density  variations  are  only 
introduced  into  the  body  force  terms  in  the  momentum  equations 

(iii)  negligible  viscous  dissipation  in  the  energy  equation 

(iv)  a  linear  dens i ty- temperature  relationship  with  a  coefficient  of 
thermal  expansion  3 

(v)  negligible  coriolis  and  centrifugal  forces  in  the  <j>  and  r  momen¬ 
tum  equations  respectively. 


The  boundary  conditions  are  as  follows.  At  the  periphery  and  the 
closed  end  of  the  cylinder  a  thermal  resistance  is  specified  based  on  the 
wall  composite  thermal  resistivity  and  an  external  convective  resistance 
based  on  empirical  correlation  of  heat  transfer  from  horizontal  cylinders 
(see  McAdams  [2]).  An  arbitrary  ambient  temperature  was  chosen.  Since  the 
open  end  is  exposed  to  hot  gases  which  will  also  be  under  the  influence  of 
convective  heat  transfer,  various  boundary  conditions  were  explored. 


5 


The  conditions  chosen  were: 


0)  a  constant  uniform  temperature  and  zero  velocity  to  represent  a 
uniform  static  environment 

( i i )  a  linear  temperature  profile  in  the  vertical  direction  and  zero 
velocity  to  represent  a  stratified,  static  environment 

(iii)  a  linear  temperature  profile  as  in  (_i  i )  but  with  a  constant 
finite  velocity  to  represent  a  convecting  downward  flow. 

In  case  (i )  the  constant  temperature  was  chosen  as  the  initial  hot  tempera¬ 
ture  of  the  gas.  The  linear  temperature  profiles  of  cases  (ii)  and  (iii) 
were  ±20°F  about  the  assumed  initial  starting  temperature.  The  relative 
effects  of  each  of  these  conditions  will  be  discussed  later. 

The  boundary  conditions,  expressed  in  terms  of  the  nondimensional 
variables  defined  in  the  Nomenclature  are: 

R"l:  WR 


0 

aw 

3R 


“Z 


30 

3R 


c)  U 

3R 


0-1 


3v 

w<i>  "  ‘  3Z 

wz  =  0 


36  = 

3z 


o2 
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Z  =  L: 


CD 


311 

“R  "  3Z 


3V 

%  _  ”  3Z 

03z  =  0  (STAGNANT-STRATI  FI  ED) 

9  =  Fl  (R,4>,Z) 

(2)  u  =  F2(R,<J>,Z)  (CONVECTIVE-STRATI  FI  ED) 

CONDITION 

e  =  F  3  (R ,  <J> ,  Z ) 

where  the  functions  Fj,  and  F3  represents  the  assumed  non dimensional  tem¬ 
perature  profiles,  and  F2  is  the  curl  of  the  assumed  velocity  profile  at 
Z  =  L. 

Mumer i cal  Procedure 

The  numerical  modelling  scheme  used  to  solve  the  system  of  Equations 
(1)  -  (5)  was  a  modified  version  of  that  used  by  Aziz  and  Heliums  [1]  for 
a  rectangular  enclosure.  The  grid  arrangement  is  shown  in  Figure  3-a.  In 
order  to  evaluate  the  dependent  variables  at  the  centerline,  a  Cartesian 
coordinate  system  was  imposed  along  the  line  r  =  0  as  indicated  in  Figure 
3-b.  This  eliminates  the  differencing  equations  from  blowing  up  along  this 
line  and  was  found  to  be  more  accurate  than  using  the  limiting  form  of  the 
differential  equations  at  r  =  0  (see  Kee  and  McKillop  [3l). 

The  numerical  procedure  consists  of  three  stages  per  time  step.  First 
the  parabolic  equations  are  solved  for  A^  which  represents  a  first 
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approximation  of  the  dependent  variable  due  to  changes  in  the  <j>  direction 

(2) 

at  a  given  time  step,  then  a  second  approximation,  A  ,  at  the  same  time 
step  due  to  changes  in  the  r  direction  is  obtained,  and  lastly  is 

calculated  based  on  changes  in  the  2  direction.  The  system  of  equations 
solved  is: 


0) 

(°). 

An+l 

A 

n 

At 

(2) 

(°) 

A  ,  " 
n+1 

A 

n 

At 

(3) 

(°) 

Vl  - 

A 

n 

Ax 

0)  (°)  (°)  (°) 

c  [H.  6,  (A  .  .  +  A  )  +  6  (A  )  +  6  (A  )  +  B] 
<p  n+ 1  n  P,  n  L  n 


0)  (°) 


(2)  (°) 


C  &  VAn+l  +  V  +  *  VAn+l  +  V  + 


6z(An)  +  B] 


(1)  (0)  (2)  (0) 
-c  &  VVi  +  A„> + 

t3>  (°) 

13  Wl  *  A„>  *  B! 


In  the  above  expressions  C  represents  a  constant  and  equals  one  when  solving 
for  the  vorticity  components  and  l/Pr  in  evaluating  the  temperatures.  The 
delta  functions,  5,  represent  the  following  combinations  of  first  and 
second  derivatives: 


„  J_  r  _  u 

<J>  R  3$  2  R  3 


,  _ii_.fu.L_i 

5R  "  3R?  R!  3R 


’  Sr  '  w  It 


The  conventional  symmetric  finite  difference  forms  for  the  first  and 
second  derivatives  were  used  to  evaluate  each  <5  to  assure  errors  of  the 
order  of  the  grid  size  squared.  As  an  illustration  of  how  Equations  (6) 
were  applied  consider  the  following  formulation  to  evaluate  the  first 
approximation  of  the  <j>  component  of  the  vorticity  at  an  interior  node 


-  %  < 


(*)  <M 

w,<i+l,j,k)  -  2<o,  (i  ,j,k)+u>(i-l  ,j,k) 


0)  (°) 

W,  0  >  j  « k)  ’  w  (ij,k) 


ufij.k)  ^(m»j»k) : 

R  2A<)>  n+1 

(°)  (°)  (°) 

(i+,’->>k)  -  2V’i’k}  +  w«i>  <M’j*k)  u(i.i.k) 

( - - R - 


cO  (i  +  1  ,j  ,k)  -  rt.  ( i  -  1  ,  j  ,  k) 


(°)  (°)  (°) 

WP  (i.j  +  l,k)  -  2w  ( i ,  j  ,k)  +  u)p  (i,j,-l,k) 


RV(i.j,k)-l 

R 


wR  (i ,j ,  +  l »k)  -wR  ('  ,j'l ,k) 


(°)  (°)  (°) 

u2  ( ' » j  » k+ 1 )  -  2uz  ( i , j ,k)  +  0  » J  *k" 1 ) 


-  W(i,j,k)< 


(°)  (°) 

“z  ( i » j  » k+ 1 )  -  i»z  (i , j  ,k-l)  • 
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(8) 


r  „  30  ,  .  2  au,R  -4 

+  I-C  ,2  “*♦  *  ST—  -V 


R  3  4> 


3U  ,  ^  3U  ,  V^’j’k^  WR 
WR  3R  Z  3R  R  i 


The  term  in  square  brackets  represents  the  B  term  in  Equations  (6).  Similar 

equations  result  for  to^  and  9,  (the  dimensionless  temperature).  Once 

(T) 

the  above  equation  is  solved  for  w  it  is  stored  and  the  equation  for 
(2)  .  .  (1)  n+1 

w'  (i,j,k)  is  solved  using  u'  in  the  8.  operator  as  indicated  in 

*  „  ^n+1  * 

n+1 

Equations  (8). 

It  is  desirable  to  keep  the  error  of  the  order  of  the  grid  size  squared. 
As  such  the  following  approximations  for  the  derivatives  at  the  boundary  were 
used : 

(|£)  *  zJm]  3A(i’j’k)  '  *>A(i,j-l,k)  +  A(!,J-2,k)] 

R=  I  *• 


(ffij)  (  2A(i,j,k)  +  5A(i,j-l,k)  -  i»A(i,j-2,k)  + 

A(i,j-3,k)J  (?) 

with  similar  expressions  for  the  derivatives  in  the  Z  direction  at  the 
closed  end. 

In  order  to  assure  the  least  error  introduced  in  the  finite  difference 
approximat ion  of  the  time  derivatives  it  is  necessary  that  B  and  the  velo¬ 
city  components  be  evaluated  at  n  +  h  or  n  +  1 .  This  requires  a  two  stage 
iteration  scheme.  First  the  vorticity  and  temperature  are  evaluated  using 
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the  old  values  of  B  from  which  the  vector  potential  and  velocity  components 
are  calculated  (this  procedure  is  outlined  further  on  and  requires  an 
iteration  procedure) .  These  updated  values  are  then  reinserted  into  the 
parabolic  set  of  equations  to  reevaluate  w  and  0  at  the  same  time  step. 

Mew  values,  of  ¥  and  V  are  then  calculated  and  again  used  to  update  to  and  9. 
This  is  continued  until  a  predetermined  convergence  criteria  is  satisfied 
(that  is,  the  fractional  change  in  the  updated  and  old  values  is  less  than 
some  tolerance).  Once  completed  the  entire  process  is  repeated  at  the 
next  time  step. 

There  is  the  further  requirement  to  specify  the  vorticity  at  the  solid 
boundaries,  R=l  and  Z=0.  These  are  obtained  indirectly  by  examining  the 
boundary  conditions  for  the  vector  potential.  Morean  (see  Reference  [l]) 
concludes  that  the  normal  derivatives  of  the  normal  V  component  at  a  solid 
surface  is  identically  zero.  Further,  the  tangential  component  of  ¥  to 
the  surface  also  must  vanish  to  satisfy  the  no-slip  conditions.  Transform¬ 
ing  these  conditions  into  the  vorticity  components  yields  the  following: 

R  =  1 : 


Z  =  0: 


aW  I 

"♦'fiT-JTW  (H'(i  k-'-k)  - 

"r  '  0 

pjll  „  ] 

u>2  -  -  ~  =  (Wj-l.k)  -  u(i  k-2 , k) ) 

(10) 

1 

"’<()  Iz  =  2 (aZ)  “  V^'  »j»k+2)) 

■ill  .1 

=  Yr  =  2HzT  ^tU^'’j,k+l)  -  U  ( i ,  j  ,  k+2 ) ) 

~  0 
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we** 


where  the  finite  difference  approximat ions  are  second-order  in  the  space 
coordinates.  This  formulation  implies,  that  new,  or  updated,  values  of  the 
velocity  components  are  required  to  evaluate  updated  vorticity  components 
at  the  boundaries. 

The  solution  scheme  for  the  elliptic  equations  for  V  is  as  follows.  A 
successive  overrelaxation  method  is  used  to  solve  the  finite  difference 
approximation  which  can  be  expressed  in  the  following  form: 


(m) 

S„  (i , j ,k)  = 


(m-1)  (m) 

“opt  rSn  (i.j.k+D  +  Sn  (i,.j,k-l) 

=  __  [ - m -  + 

P 


(m-1) 


(m) 


sn  (i,j+i,k>  +  sn  (i , j- 1 ,k) 


TP 

(m) 

(m) 

Sn  <5 

,j+l,k)  -  S  (i,j-l,k) 

2RAR 

(m-1 ) 

(m) 

s 

n 

(i+i,j,k)  +  sn  (i- 1 , j ,k) 

(m-  I ) 


(m-1) 


4P  Sn  -“optS  +  Sn  (T,j,k)  (II) 


where  A^  is  dependent  on  the  grid  size  and  is  defined  in  the  nomenclature. 

The  superscripts  (m)  and  (m-1)  represent  current  and  previous  iteration 
values  respectively.  The  term  §.  is  the  appropriate  component  of  vorticity 
to  correspond  with  the  particular  vector  potential  component  being  evaluated. 
In  this  expression  u>q  represents  an  acceleration  parameter  and  can  take 
on  values  typically  in  the  range  zero  to  two.  An  optimal  value  was  found 
which  is  designed  to  provide  the  fastest  rate  of  convergence.  (See  Forsythe 


and  Wason  [4] . ) 


Equation  (10)  is  iterated  upon  until  the  fractional  change  in  the  vec¬ 
tor  potential  component  is  less  than  some  set  tolerance.  The  convergence 
is  tested  at  each  nodal  location  in  order  to  assure  that  the  entire  system 
is  well  behaved.  This  was  later  changed  to  a  selected  representative  point 
to  shorten  the  computational  time. 

The  grid  spacing  was  9  x  9  x  11.  That  is,  the  radial  increment  equals 
rQ/8,  the  angle  increment  tt/8  and  the  longitudinal  increment  L/(10). 

Later  a  9  x  9  x  6  grid  system  proved  to  be  effective. 

The  program  was  run  until  steady  state  conditions  were  obtained.  Steady 
state  is  determined  by  evaluating  the  total  heat  loss  rate  from  the  cylinder. 
When  the  heat  loss  rate  no  longer  varies  within  one  percent, steady  state  is 
assumed. 

A  flow  chart  for  the  solution  algorithm  is  shown  in  Figure  L\. 

Results 

The  numerical  experiments  run  to  date  are  here  summarized.  The  solid 
boundary  thermal  resistance  was  calculated  using  a  three  layer  composite 
wall  consisting  of  refractory  brick,  an  insulating  wall  and  a  pressure 
vessel.  Values  of  the  thermal  properties  were  obtained  from  an  existing 
test  facility  at  Arnold  Air  Station,  TN.  Variation  of  these  values  by  ±20? 
did  not  alter  the  results  presented  herein.  The  relatively  low  conduc¬ 
tivity  wall  of  the  system  results  in  rather  quickly  established  steady- 
state.  However,  the  wall  was  assumed  to  be  in  local  equilibrium  with  the 
instantaneous  gas  state  which  would  tend  to  reduce  the  transient  times 


somewhat.  This  effect  was  not  rigorously  studied  in  this  investigation, 
but  an  examination  of  the  wall  thermal  capacitance  shows  that  it  should 
be  a  minor  effect. 

The  majority  of  the  results  are  presented  in  terms  of  the  nondimen- 
sional  convective  film  coefficient  -  the  Nusselts  number,  Nu  defined  as 


T  =  average  wall  temperature. 

w 


Local  values  of  Nu  were  also  found  at  each  point  on  the  wall.  Average 
values  at  each  axial  location  were  calculated  based  on  the  radial  plane 
average-  bulk  temperature  and  wall  temperature.  The  wall  heat  flux  was 
determined  at  each  point  from  the  local  wall  temperature  gradient  using  a 
three  point  approximat ion  for  the  finite  difference  representation. 

The  transient  nature  of  the  wall  heat  flux  is  best  shown  in  the  Nu  vs. 
nondimensional  time  plot  of  Figure  5-  These  results  are  for  the  open  end 
boundary  condition  case  (i)  but  all  three  cases  essentially  yield  the  same 
results.  Steady  state  conditions  occur  when  t  ^  .005.  The  experimental 
results  of  Evans  and  Stefany  [5]  for  transient  convection  in  a  horizontal 
cylinder  with  a  step  change  in  wall  temperature  indicate  steady  values  of 
the  film  coefficient,  h,  after  approximately  20  seconds  using  n-butanal . 
Evaluating  the  thermal  properties  of  n-butanal  at  the  average  fluid  tempera¬ 
ture  used  by  Evans  and  Stefany  in  their  2.5  inch  diameter  cylinder  results 
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in  a  value  of  t  =  .00^8  which  is  in  agreement  with  these  numerical  results. 

The  boundary  condition  imposed,  of  uniform  wall  cooling  results  in 
cooler  fluid  near  the  boundary  than  in  the  central  region.  The  cooler  fluid 
sinks  and  sets  up  a  clockwise  convection  current.  This  motion  is  depicted 
in  Figure  6.  Typical  isotherms  and  two  dimensional  vector  potentials 
(stream  lines)  are  shown  as  functions  of  time.  An  upwardly  buoyed  motion 
is  seen  by  the  mushroomed  shaped  isotherms.  The  vortex  type  motion  has  a 
sinking  center  of  rotation  with  increasing  time. 

The  induced  motion  results  in  only  slight  variations  of  the  local  Mu 
around  the  cylinder.  Figure  7  shows  the  increased  heat  transfer  rates  near 
the  top  of  the  cylinder. 

Axial  variations  in  local  values  of  Nu  are  shown  in  Figure  8  for  all 
three  cases  studied.  At  Z=0  the  solid  boundary  impedes  the  convective 
motion  resulting  in  lower  values  of  Hu.  At  the  open  boundary,  Z=L,  only 
slightly  higher  values  of  Hu  result  for  case  (iii)  compared  to  (i)  and  (ii) 
where  case  (iii)  is  the  imposed  constant  velocity  and  linear  temperature 
gradient.  Removing  the  velocity  but  retaining  the  temperature  gradient 
also  slightly  increases  the  radial  heat  transfer  from  the  cylinder.  In 
addition,  the  axial  location  of  maximum  heat  transfer  is  shifted  back 
towards  the  closed  end.  A  rough  schematic  of  the  axial  isotherm  and  vector 
potential  lines  are  shown  in  Figure  9  for  case  (i).  A  hot  central  core 
exists  with  temperature  gradients  toward  the  conducting  surfaces.  Two  longi 
tudinal  convection  cells  exist  (for  the  range  of  Ra  investigated). 

A  correlation  of  the  Hu  vs.  Ra  was  calculated  for  the  geometric  and 
boundary  condition  constraints  imposed.  The  Rayleigh  number  is  defined  in 
terms  of  the  temperature  difference  (T,  -  T  )  and  fluid  properties  were 
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figure  8.  Longitudinal  variation  of  the  local  Nusselt  number  for 
the  three  boundary  conditions  imposed  at  the  open  end. 
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evaluated  at  T  .  Correlations  of  the  type  Nu  =  CRn  are  presented  in  Table  1 

D  3 

for  various  geometries  and  boundary  conditions  found  in  the  literature. 

These  type  of  correlations  tend  to  smooth  out  the  boundary  conditions 
imposed  as  well  as  the  system  geometry.  The  present  study  results  are 
given  in  the  table  for  n  =  .25  where  the  length  scale  is  chosen  to  be  the 
cylinder  diameter  (entree  6).  For  comparison  with  Haas  [8],  the  uniform 
heat  flux  case,  a  .21  power  was  forced  with  the  resulting  value  of  C  =  1.81 
which  is  kS%  higher  than  that  found  by  Maas.  Changing  the  length  scale,  L, 
in  evaluating  Nu  and  Ra  to  represent  three  dimensional  effects  such  that 
l  m  1_  J_ 

L  LH  ld 

where  L„  is  the  horizontal  axial  length  and  L_  the  diameter  yields  the  result 
of  entree  8  in  the  table.  Although  this  is  closer  to  the  correlations  of 
Maas  and  Deaver  and  Eckert  higher  values  of  Nu,  for  the  same  Ra  are  calculated. 
This  is  suspected  to  be  a  result  of  the  end  boundary  effects  -  the  additional 
induced  buoyancy  in  the  flow  which  is  shown  in  Figure  9  may  enhance  the  heat 
transfer  rates,  compared  to  the  two  dimensional  case. 

Representative  values  for  the  natural  convective  film  coefficient  were 
incorporated  into  an  existing  finite  element  radiation  model  to  determine  wall 
temperature  values.  This  was  carried  out  at  the  Arnold  Air  Station,  TN  and 
applied  to  a  specific  test  facility.  To  date,  experimental  verification  of 
the  results  have  not  been  complete.  The  analysis  indicates  an  approximate 
ten  percent  increase  in  the  overall  wall  heat  transfer  rate  due  to  natural 


convection. 


TABLE  1.  Comparison  of  Mean  Nusselt  Number  Correlations 


Investigator* 

n 

£ 

Cond i t ion 

1. 

Deaver  £  Eckert  [7] 

.214 

1.181 

Experiment,  increasing 
wall  temperature 

2. 

Maas  [8] 

.210 

1.215 

Experimental,  uniform 
wal 1  heat  flux 

3. 

Evans  £  Stefany  [9] 

•  25 

0.55 

Experimental,  step  change 
in  wall  temperature 

4. 

Kuehn  £  Goldstein  [10] 

LT\ 

CM 

0.20 

Experimental,  annulus 

5. 

Ozoe,  Sayama  and 

Church  ill  [  1 1  ] 

•  336 

.0981 

Experimental,  rectangular 
channel,  aspect  ratio  =  1 

6. 

Present  Study 

.25 

1.13 

Numerical,  three  dimensional 

7. 

Present  Study 

,210 

1.81 

Numerical,  three  dimensional 

8. 

Present  Study 

.210 

1.61 

Lenqth  scale  =  T~  +  r~ 

H  LR 

•'Sources  are  listed  in  Reference  section 
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RECOMMENDATIONS  FOR  FURTHER  STUDIES 


The  numerical  model  developed  in  this  project  was  applied  to  a  specific 
configuration  (a  fixed  length-to-diameter  ratio)  and  boundary  conditions 
(constant  wall  heat  flux).  Only  a  limited  range  of  Ra  was  investigated  to 
be  consistent  with  the  Arnold  Air  Station  test  facility  conditions.  Only 
slight  variations  in  axial  heat  flux  was  found,  compared  to  two  dimensional 
mode  Is. 

If  larger  values  of  Ra  are  to  be  imposed,  which  may  be  the  case  for 
very  high  gas  temperatures  or  if  using  a  high  Prandtl  number  fluid  a  transi¬ 
tion  to  turbulent  conditions  will  result.  It  is  recommended  that  in  order 
to  study  this  problem  the  numerical  model  be  generated  in  two  dimensions 
(r  and  <(>)  so  that  a  much  finer  grid  space  can  be  used.  Turbulence  modelling 
should  then  be  introduced,  such  as  a  simple  algebraic  stress  formulation  or 
a  more  rigorous  turbulent  transport  analysis  such  as  has  been  applied  to 
external  buoyant  shear  flows  by  Liburdy,  Groff  and  Faeth  [5]  and  Liburdy 
and  Faeth  [6], 

In  addition,  in  order  to  test  these  models,  it  is  necessary  to  carry 
out  experimental  verifications.  These  studies  should  be  capable  of  mea¬ 
suring  surface  heat  flux  rates,  bulk  fluid  temperatures,  and  velocity  and 
temperature  profiles.  In  addition,  for  extension  into  the  turbulent  regime, 
fluctuating  temperature  and  velocity  components,  and  appropriate  cross 
correlations  as  possible  should  be  measured.  Only  a  limited  number  of  such 
studies  have  been  performed  in  natural  convecting  flow  in  enclosures  in 
general.  The  instrumentation  must  be  rather  sophisticated  in  order  not  to 
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interrupt  the  flow  field  (laser  velocimeter  techniques  seem  most  appropriate 
for  velocity  measurements  and  fine  wire  thermocouples  for  temperature  data 
have  been  used  previously  to  yield  satisfactory  resultsl. 
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PUBLICATIONS  AND  PAPER  PRESENTATIONS 

Some  of  the  information  contained  in  this  report  which  has  resulted 
from  U.S.A*F.  Contract  F49620-79-C-01 68  has  been  presented  at  the  South¬ 
eastern  Thermal  Science  Seminar,  Florida  Atlantic  University,  Boca  Raton, 
Florida,  May,  1979  and  appears  in  the  Abstracts  of  that  meeting. 

Further  modifications  to  the  computational  scheme  in  future  studies 
as  well  as  matching  experimental  data  is  anticipated  as  discussed  in  the 
section  on  Recommendations  for  Further  Study.  Such  future  investigations 
are  expected  to  add  to  the  rather  limited  archival  knowledge  of  turbulent 
natural  convection  in  enclosures. 
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